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Abstract 

In this article we derive an unbiased expression for the expected mean-squared error 
associated with continuously differentiable estimators of the noncentrality parameter of a chi- 
square random variable. We then consider the task of denoising squared-magnitude magnetic 
resonance image data, which are well modeled as independent noncentral chi-square random 
variables on two degrees of freedom. We consider two broad classes of linearly parameterized 
shrinkage estimators that can be optimized using our risk estimate, one in the general context 
of undecimated filterbank transforms, and another in the specific case of the unnormalized Haar 
wavelet transform. The resultant algorithms are computationally tractable and improve upon 
state-of-the-art methods for both simulated and actual magnetic resonance image data. 
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I. Introduction 

Magnetic resonance (MR) imaging is a fundamental in vivo medical imaging technique that 
provides high-contrast images of soft tissue without the use of ionization radiation. The signal-to- 
noise ratio (SNR) of an acquired MR image is determined by numerous physical and structural 
factors, such as static field strength, resolution, receiver bandwidth, and the number of signal 
averages collected at each encoding step [1], [2]. Current developments in MR imaging are 
focused primarily on lowering its inherently high scanning time, increasing the spatiotemporal 
resolution of the images themselves, and reducing the cost of the overall system. However, the 
pursuit of any of these objectives has a negative impact on the SNR of the acquired image. A 
post-acquisition denoising step is therefore essential to clinician visualization and meaningful 
computer-aided diagnoses [1]. 

In MR image acquisition, the data consist of discrete Fourier samples, usually referred to as 
/c-space samples. These data are corrupted by random noise, due primarily to thermal fluctuations 
generated by a patient's body in the imager's receiver coils [1], [2]. This random degradation 
is well modeled by additive white Gaussian noise (AWGN) that independently corrupts the real 
and imaginary parts of the complex-valued fc-space samples [3]. Assuming a Cartesian sampling 
pattern, the output image is straightforwardly obtained by computing the inverse discrete Fourier 
transform of the fc-space samples; the resolution of the reconstructed image is then determined by 
the maximum fc-space sampling frequency. Note that some recent works (see [4] and references 
therein) aim to accelerate MR image acquisition time by undersampling the /c-space. Non- 
Cartesian (e.g., spiral or random) sampling trajectories, as well as nonlinear (usually sparsity- 
driven) reconstruction schemes are then considered. This axis of research is, however, outside 
the scope of the present work. 

Following application of an inverse discrete Fourier transform, the resulting image may be 
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considered as a complex-valued signal corrupted by independent and identically distributed 
samples of complex AWGN. In magnitude MR imaging, the image phase is disregarded and 
only the magnitude is considered for visualization and further analysis. Although the samples of 
the magnitude image remain statistically independent, they are no longer Gaussian, but rather are 
Rician distributed [5]. Contrary to the Gaussian case, this form of "noise" is signal-dependent 
in that both the mean and the variance of the magnitude samples depend on the underlying 
noise-free magnitude image. Consequently, generic denoising algorithms designed for AWGN 
reduction usually do not give satisfying results on Rician image data. 

Denoising of magnitude MR images has thus gained much attention over the past several 
years. Two main strategies can be distinguished as follows. In the first, the Rician data are 
treated directly, often in the image domain. In the second, the denoising is applied to the squared 
magnitude MR image, which follows a (scaled) noncentral chi-square distribution on two degrees 
of freedom, whose noncentrality parameter is proportional to the underlying noise-free squared 
magnitude. An appealing aspect of this strategy is that it renders the bias due to the noise 
constant rather than signal dependent, thus enabling standard approaches and transform-domain 
shrinkage estimators to be applied. 

Following the first strategy previously mentioned, several Rician-based maximum likelihood 
estimators [6]-[8] have been proposed. In [8], this estimation is performed non-locally among 
pixels having a similar neighborhood. A Bayesian maximum a posteriori estimator has been 
devised in [9], with the prior modeled in a nonparametric Markov random field framework. 
Very recently, Foi proposed a Rician-based variance stabilizing transform which makes the use of 
standard AWGN denoisers more effective [10]. Following the second strategy, a linear minimum 
mean-squared error filter applied to the squared magnitude image has been derived in [11], [12]. 
In addition to these statistical model-based approaches, much work has also been devoted to the 
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adaptation and enhancement of the nonlocal means filter originally developed by Buades etal. 
for AWGN reduction [13]. The core of this relatively simple, yet effective, denoising approach 
consists of a weighted averaging of similarly close (in spatial and photometric distance) pixels. 
Some of these adaptations operate directly on the Rician magnitude image data [14], [15], while 
others are applied to the squared magnitude image [16]— [18]. 

In addition to these image-domain approaches, several magnitude MR image denoising algo- 
rithms have also been developed for the wavelet domain. The sparsifying and decorrelating effects 
of wavelets and other related transforms typically result in the concentration of relevant image 
features into a few significant wavelet coefficients. Simple thresholding rules based on coefficient 
magnitude then provide an effective means of reducing the noise level while preserving sharp 
edges in the image. In the earliest uses of the wavelet transform for MR image denoising [19], 
[20], the Rician distribution of the data was not explicitly taken into account. Nowak subsequently 
proposed wavelet coefficient thresholding based on the observation that the empirical wavelet 
coefficients of the squared-magnitude data are unbiased estimators of the coefficients of the 
underlying squared-magnitude image, and that the residual scaling coefficients exhibit a signal- 
independent bias that is easily removable [21]. While the pointwise coefficient thresholding 
proposed in [21] is most natural in the context of an orthogonal discrete wavelet transform, 
Pizurica etal. subsequently developed a Bayesian wavelet thresholding algorithm applied in an 
undecimated wavelet representation [22]. Wavelet-based denoising algorithms that require the 
entirety of the complex MR image data [23]-[25] are typically applied separately to the real and 
imaginary components of the image. 

In this article, we develop a general result for chi-square unbiased risk estimation, which we 
then apply to the task of denoising squared-magnitude MR images. We provide two instances 
of effective transform-domain algorithms, each based on the concept of linear expansion of 
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thresholds (LET) introduced in [26], [27]. The first class of proposed algorithms consists of a 
pointwise continuously differentiable thresholding applied to the coefficients of an undecimated 
filterbank transform. The second class takes advantage of the conservation of the chi-square 
statistics across the lowpass channel of the unnormalized Haar wavelet transform. Owing to 
the remarkable property of this orthogonal transform, it is possible to derive independent risk 
estimates in each wavelet subband, allowing for a very fast denoising procedure. These estimates 
are then used to optimize the parameters of subband-dependent joint inter-/intra-scale LET. 

The article is organized as follows. We first derive in Section II an unbiased expression 
of the risk associated with estimators of the noncentrality parameter of a chi-square random 
variable having arbitrary (known) degrees of freedom. Then, in Section III, we apply this result 
to optimize pointwise estimators for undecimated filterbank transform coefficients, using linear 
expansion of thresholds. In Section IV, we give an expression for chi-square unbiased risk 
estimation directly in the unnormalized Haar wavelet domain, and propose a more sophisticated 
joint inter-/intra-scale LET. We conclude in Section V with denoising experiments conducted 
on simulated and actual magnitude MR images, and evaluations of our methods relative to the 
current state of the art. Note that a subset of this work (mainly part of Section III) has been 
accepted for presentation at the 2011 IEEE International Conference on Image Processing [28]. 

II. A Chi-Square Unbiased Risk Estimate (CURE) 

Assume the observation of a vector y £ of N independent samples y n , each randomly 
drawn from a noncentral chi-square distribution with (unknown) noncentrality parameter x n > 
and (known) common degrees of freedom K > 0. We use the vector notation y ~ x|-( x )> 
recalling that (for integer K) the joint distribution p(y|x) can be seen as resulting from the 
addition of K independent vectors on M.+ whose coordinates are the squares of non-centered 
Gaussian random variables of unit variance. This observation model is statistically characterized 
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N N , x 

P{y\^) = X\p(Vn\x n ) = \\^"' M ^ [—) I*-i(y/x n y n ), (1) 
n=l n=l 

El /^\2fc+a 
— — - - is the a-order modified Bessel function of the 

£;!r(£: + a + 1) \2J 

first kind. 

The chi-square distribution of (1) is most easily understood through its characteristic function, 
or Fourier transform: 



,rX r . 

N exp 



) 

For instance, by equating first and second order Taylor developments of this Fourier transform 
with the corresponding moments, it is straightforward to show that 

#{y} = x + K-l, 
and g { ||y || 2 } = ||x|| 2 + 2(K + 2)l T x + NK(K + 2), 
where $ {■} is the expectation operator. 

When designing an estimator x = f (y ) of x, a natural criterion is the value of its associated 
risk or expected mean-squared error (MSE), defined here by 

f 1 i 1 N 

S {MSE} = S\- ||f (y) - x|| 2 = - {& {/"(y) 2 } " 2^ {xnUy)} + x 2 n ) . (4) 

^ > n=l 

The expectation in (4) is with respect to the data y; the vector x of noncentrality parameters 
may either be considered as deterministic and unknown, or as random and independent of y. 

In practice, for any given realization of the data y, the true MSE cannot be computed because 
x is unknown. Yet, paralleling the general case of distributions in the exponential family [29], 
it is possible to establish a lemma that allows us to circumvent this issue, and estimate the risk 
without knowledge of x. The main technical requirement is that each component / n (y) of the 
estimator f : M N — > R N be continuously differentiable with respect to y n , and so we introduce 
the following notation: 0f(y) = [j-/„(y)l , <9 2 f(y) = \&f n (y) 



Kn<N 
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Lemma 1. Assume that the estimator f(y) = [/ n (y)]i<n</v is such that each f n (y) is continu- 
ously dijferentiable with respect to y n , with weakly differentiable partial derivatives df n (y ) / dy n 
that do not increase "too quickly" for large values of y; i.e., there exists a constant s < 1/2 
such that for every 1 < n < N, lim^^+oo e~ syn | g|-/n(y)| = 0. Then 

S {x T f (y)} = & {(y — K ■ l) T f (y)} - AS j (y - | ■ l)' 9f (y) | + 4 <T {y T a 2 f (y)} . (5) 

Proof: We first evaluate the scalar expectation <f {x„/ n (y)} appearing in (4), before sum- 
ming up the contributions over n to get the expectation of the inner product x T f(y). 

Let us consider the characteristic function given by (2). By differentiating the logarithm of 
this Fourier transform w.r.t. the variable to n , we find that 

1 dp(u\x) _ jK jx n 
p(u>|x) duj n l + 2juj n (l + 2jw n ) 2 ' 

After a rearrangement of the different terms involved, we get 

x n p(u\x) = j(l + 2ju n f dp ^ ) - K(l + 2jw n )p(w|x). 

OiO n 

Using (1 + 2juj n ) 2 = 1 + 4juj n + 4(ju} n ) 2 and recalling that multiplication by joj n corresponds 
to differentiating w.r.t. y n , while differentiation w.r.t. uj n is equivalent to a multiplication with 
— jy n , we can deduce that the probability density p(y|x) satisfies (in the sense of distributions) 
the following linear differential equation 

won.) = (i + 4A + 4|i){ ! ,„ P ( y |x,} - K(i +2 A) W yi*)}. 

The expectation $ {x n fn{y)} is simply the Euclidean inner product (with respect to Lebesgue 
measure) between / ra (y) and x„p(y|x). Continuous differentiability of f(y) coupled with weak 
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differentiability of <9f(y) implies that we may twice apply integration by parts to obtain 

& {Xnfn{y)} = (fn(y),X n p(y\x)) 

= <>"&* i 1 + 'Wn + 2 4 !g) {! "' P<y|X)} ) - A '( ; " (y) ' + <P<y|X)} 



2 / dy n dyl 
plus additional integrated terms that do not depend on x n , which vanish because of our (con- 
servative) assumption on how / n (y) increases when y n — > +00. (Asymptotic expansions of 
Bessel functions can be used to show that whenever s < 1/2, we have that e sy "p(y|x) — > 
and e syn dp(y\x) /dy n — > 0.) Finally, summing up over the index n yields (5). ■ 
We can then deduce a theorem that provides an unbiased estimate of the expected MSE given 
in (4). We term this random variable CURE, an acronym for Chi-square Unbiased Risk Estimate. 

Theorem 1. Let y ~ Xk( x ) an d assume that f (y) satisfies the regularity conditions of Lemma 1. 
Then, the following random variable: 

1 ( »- - > "-"i2 .1/.. k 1 ^^ , 8 ((,. k i\ T s*u.\ -T»2- 



cure = - ^|f( y )-( y -^.i)^-4 T (y-yiJJ + jj {{y—2' 1 ) 9f (y)-y T<5 f (y)J- ( 6 ) 

is an unbiased estimate of the risk; i.e., £ {CURE} = $ {MSE}. 

Proof: As with other unbiased risk estimates, we express the MSE as a sum of three terms: 
||f(y) - x|| 2 = \\mf- 2£f(y) + \\xf , 

term 1 term 2 term 3 

which we replace by a statistical equivalent that does not depend on x anymore. 

Term 1 needs no change, term 2 can be replaced according to (5) from Lemma 1, and term 3 
can be reformulated using the noncentral chi-square moments of (3) to yield 

IWI 2 = 8 { lly II 2 } " 2(K + 2)1 T ^ {y} + NK(K + 2). 
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Putting everything together then leads directly to (6), and thus the theorem is proved. ■ 
In Sections III and IV we propose specific examples of CURE-optimized transform-domain 
processing for estimating the unknown noncentrality parameter vector x from data y. 

III. CURE-Optimized Denoising in Undecimated Filterbank Transforms 

In this section, we focus on processing within a ( J + l)-band undecimated filterbank trans- 
form, as depicted in Fig. 1. This broad class of redundant representations notably includes the 
undecimated wavelet transform (UWT) and overlapping block discrete cosine transform (BDCT). 
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Fig. 1. Undecimated ( J + l)-band analysis/synthesis filterbank. 



A. Image-domain CURE for transform-domain processing 

Nonlinear processing via undecimated filterbank transforms is a general denoising strategy 
long proven to be effective for reducing various types of noise degradations [27], [30]-[34]. It 
essentially boils down to performing a linear (and possibly redundant) analysis transformation 
of the data, which provides empirical coefficients that are then thresholded (possibly using a 
multivariate nonlinear function), the result of which is finally passed to a linear synthesis trans- 
formation. When treating signal-dependent noise, the entire denoising procedure is conveniently 
expressed in the generic form f(y) = R0(Dy, Dy) [34], where 
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• The circulant matrices D = [Dq D{ . . . Dj] T and R = [Ro Ri . . . Rj] implement an 
arbitrary pair of analysis/synthesis undecimated filterbanks (Fig. 1) such that the perfect 
reconstruction condition RD = Id is satisfied. Each component of the N x N circulant 
submatrix = [d{j]i< k ,i<N and R^ = [r J k;l ]i<k,i<N is given by 
' d i,i = T,9^-k + nN] 



, forj = 0... J. (7) 

i,l = 229j[k-l + nN] 

We assume that the considered analysis filters have unit norms; i.e., ^ n (d\ n ) 2 = 1, for 
I = 1 . . . N, j = . . . J. By convention, we also assume that, for j = 1 ... J, each Dj 
implements a highpass channel; i.e., V/, ^2 n dj n = 0. The complementary lowpass channel 
is then implemented by D with VI, ^ n d^ n = 2 J / 2 . 

Hence, denoting by Wj = D^y = [wj]i<i<N (resp. u>j = D^x = [u^Jixkat) each vector 

of noisy (resp. noise-free) transform coefficients, we have 

sivA\ = uL for ?' = 1 ... J, 
V/ = 1 . . . AT ,{ I ' J ' (8) 

While the noisy highpass coefficients are unbiased estimates of their noise-free counterparts, 
the lowpass coefficients exhibit a constant bias (2 J / 2 K) that must be removed. 

• The circulant matrix D = [Dq . . . Dj] T implements a linear estimation of the variance 
Wj = Djy = [wj]i<i<N of each transform coefficient wj. The actual variance is given by 

N N , x 

var {wj} = ^(4 n ) 2 var{y n } ( => 4 J>( n ) 2 (x n + - . (9) 

n=l n=l ^ ' 

Since £ {y n } = x n + K, the natural choice Dj = [(dj n ) 2 ]i<i,n<N achieves 

H}=4(^{^}-f). (10) 

• The vector function : R L x R L ->■ R L , where L = ( J + 1)N, can generally be arbitrary, 
from a simple pointwise thresholding rule to more sophisticated multivariate processing. In 
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this section, we will however consider only subband-adaptive pointwise processing; i.e., 



0(w,w) = [Oi(wj,wj)] <j<j,i<i<N- 



(11) 



We further assume that the transform-domain pointwise processing of (11) is (at least) contin- 
uously differentiable, with piecewise-differentiable partial derivatives. Introducing the notation 



<9i#(w, w) 



<96^(w, w) 



dwi 



, <9 2 0(w, w) 



d^O^w, w) = 



d 2 9i(w, w) 
dwf 



1<1<L 
|2 



d6i(w, w) 



dwi 



,<9 22 0(w, w) = 



J KKL 



d 2 6i(w, w) 



9W 1 J KKL 



KKL 
|2 



,a 12 0(w,w) 



d 2 0i(w, w) 
dwidwi 



and denoting by "0" the Hadamard (element-wise) matrix product, we have the following. 

Corollary 1. For pointwise processing of the form given by (11) and satisfying the requirements 
of Lemma 1, the risk estimate of (6) takes the following form: 



CURE 



1 

N 



||f(y)-(y-K-l)|| 2 -4 T 



K 



+ Jj ( y " f ' *) ( (R DT ) 9l0 ( w ' w) + (R D T )<9 2 0(w, w)) 
y T ((R D T D T )9 1 2 1 6>(w, w) + (R D T D T )d 2 2 0(w, w)) 



N 
16 



+ ^y T (R D T D T )d 2 2 0(w, w). 



(12) 



The proof of this result is straightforwardly obtained by developing the term (y—K/2 ■ l) T <9f (y)- 
y T 9 2 f(y) from (6) for 0(w,w) as defined in (11). A similar result for transform-domain 
denoising of mixed Poisson-Gaussian data is proved in [34]. 

For the remainder of this section, we drop the subband superscript j and the in-band location 
index I. We thus denote by w, w, and uj any of the wj, wj, and wj, for j = 1 . . . J. 
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B. Choice of thresholding rule 

We need to specify a particular shrinkage or thresholding rule to estimate each unknown 
highpass coefficient oj from its noisy counterpart w. In the minimum-mean-squared error sense, 
the optimal pointwise shrinkage factor is given by 

• fi f ( x2 l (8),d0) 1 4 (g {W} - f ) 

a = arg mm & \{aw - oj) \ = 1 — — . (13) 

There are various possible implementations of the above formula to yield an effective shrinkage 
function. For the case of K = 2 degrees of freedom, Nowak proposed in [21] the function 

„/ / > 4max(u; — 1, 1) \ 

9(w,w) = max (l - A ^ — ,0ji«, (14) 

where the particular choice A = 3 was motivated by a Gaussian prior on the noisy coefficients 
w. Our experiments have indicated that replacing max(W— 1, 1) by w gives slightly better MSE 
performance, and so, following the recent idea of linear expansion of thresholds (LET) [27], we 
propose the following shrinkage rule for arbitrary degrees of freedom: 

9(w,w;a) = y^q^max — Xj— ^,o)w, with I « N, (15) 

i=i N v—^- ' 

9i(w,w) 

which can be seen as an optimized generalization of (14). To satisfy the requirements of 
Corollary 1, we implement a continuously differentiable approximation of the max(-) function. 

Empirically we have observed 1 = 2 terms per subband to be the best choice in (15). The 
vector a G R 1 of subband-adaptive parameters can be optimized in closed form via least squares, 
while Ai and A2 can be optimized by minimizing the risk estimate of (12) directly. However, 
fixing Ai = 3 and A2 = 9 was observed to work well in all of our experiments, and leads to 
a much faster implementation. (We observed values close to (3,9) to yield equivalent results 
±0.2 dB.) A potential realization of the proposed LET is displayed in Fig. 2. 
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Fig. 2. Possible realization of the proposed thresholding rule of (15) (7 = 2, Ai = 3, A2 = 9, a = [0.75 0.25] T ). 



C. Implementation 

The overall subband-adaptive transform-domain estimator is thus, for Z = {0 ... J} x {1, 2}, 

J 2 

f (y) = £ ai fi(y) = a? (R w - K) + £ £ afo-^w,-, w,-), 



Lowpass 
bias removed 



j=l i=l 



with the CURE-optimized parameters a = [a^gx the solution to Ma = c, where 

K 



(y-K- 1) T f,(y) - 4 ( ( y - - • 1 ) ^(y) - y T ^fi(y) 



J iei 



(16) 



( M = [fi(y) T fj(y)] lj6l . 

IV. CURE-Optimized Denoising via Unnormalized Haar Wavelet Transform 

In the previous section, we have considered the general case of an undecimated filterbank 
transform and derived the corresponding image-domain MSE estimate. Owing to the intractability 
of the noncentral chi-square distribution after an arbitrary (even orthogonal) transformation, an 
explicit transform-domain risk estimate is generally unobtainable. Remarkably, in the particular 
case of the unnormalized Haar wavelet transform, the derivation of such an explicit subband- 
dependent MSE estimate is possible. Its construction is presented in this section. 
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Same scheme 

applied recursively i 
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Fig. 3. Signal-dependent noise reduction in the unnormalized Haar discrete wavelet transform. 



A. Unnormalized Haar wavelet-domain CURE 

The ID unnormalized Haar discrete wavelet transform consists of a critically- sampled two- 
channel filterbank (see Fig. 3). On the analysis side, the lowpass (resp. highpass) channel 
is implemented by the unnormalized Haar scaling (resp. wavelet) filter whose z-transform is 
H{z~ x ) = 1 + (resp. G(z~ l ) = 1 — z~ x ). To achieve a perfect reconstruction, the synthesis 
side is implemented by the lowpass and highpass filters H{z) = (l + z)/2 and G(z) = (l — z)/2. 
At a given scale j, the unnormalized Haar scaling and wavelet coefficients of the observed data 
y = s° are given by 



j-i i j'-i 
+ s J 



n — °2n "T °2n-l> W n ~ °2n 



'2n-l- 



(17) 



Similarly, the unnormalized Haar scaling and wavelet coefficients of the noncentrality param- 
eter vector of interest x = are given by 



j j — l _. j— l 

— ?2n + ? 2n-l 



? 2n ? 2n-l- 



(18) 



Since the sum of independent noncentral chi-square random variables is a noncentral chi- 
square random variable whose noncentrality parameter and number of degrees of freedom are 
the summed noncentrality parameters and number of degrees of freedom [35], the empirical 
scaling coefficients follow a noncentral chi-square distribution; i.e., s- 7 ~ x|- where Kj = 
2 J K. Moreover, since the squared lowpass filter coefficients are the lowpass filter coefficients 
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themselves, the scaling coefficients can be used as estimates of the variance of the (same-scale) 
wavelet coefficients. In the notation of Section III-A, this means that w J = s 3 for j = 1 ... J. 

Denoting by Nj the number of samples at a given scale j and assuming a subband-adaptive 
processing 6 3 : R N3 x R N3 ->• R Ni , the MSE in each highpass subband j is given by MSEj = 
■^-\\0 J (w 3 , s J ) — oj j || 2 , and we have the following theorem. 

Theorem 2. Let 0(w, s) = 9 3 : (w J s 3 ') be an estimator of the unnormalized Haar wavelet 
coefficients u> = u> 3 of x at scale j, satisfying the conditions of Lemma 1. Then the random 
variable 

CURE, = ±- (j|0(w, s) - w|| 2 -4 T (s- ^ • l) ) + ± ( (s- ^ • I)" ft0(w, s) + w T d 2 0(w, s)) 

- J- (w T (^(w, s) + di 2 0(w, s)) + 2s T d 2 2 0(w, s)) (19) 
is an unbiased estimate of the risk for subband j; i.e., £ {CURE, } = $ {MSE, }. 



Proof: We consider the case j = 1, so that we may use K = Kj/2, y = s 3 1 , and x = q 3 1 
to ease notation. We first develop the squared error between u and its estimate 0(w, s): 

S (||6>(w,s) -cj|| 2 } = S { ||6>(w, s)|| 2 } - 2£ {u> T 0(w,s)} + \\u\\ 2 . (20) 

(i) (ii) 

We can then evaluate the two expressions (I, II) that involve the unknown u>. 

(I) Computation of $ {u> T 6(w, s)} = ^ & {u n 9 n (w, s)}: We can successively write 

n=l 

(18) 

(f {W„(9 n (w,s)} = # {x 2 rA(w,s)} - S {x 2n _i0 n (w, s)} 

(5) ± ll) g {w n (9 n (w, s) + 4 (a 2 ^ n (w, s) + <9 2 A(w, s) - d 2 #„(w, s))) } 
-4<?{(s„ - K)di0 n (w,8) - 2s n df 2 e n (w,s)} . (21) 
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AT, 

(II) Computation of ||cj|| 2 = ^w 2 : We can successively write 

71=1 

oj 2 n = S {cj n w n } ( = ] £{w 2 n - 4(s n -K)}. (22) 

Inserting (21), (22) into (20) yields the desired equality; for j > 1, the proof is similar. ■ 
Subband superscript j will be omitted below, as we consider any of the J wavelet subbands. 

B. CUREshrink 

A natural choice of subband-adaptive estimator in orthogonal wavelet representations is soft 
thresholding, introduced by Weaver etal. [19] and theoretically justified by Donoho [36]. In 
contrast to the AWGN scenario, a signal-dependent threshold is required here. As in the case of 
Poisson noise removal [37], [38], we wish to adapt the original "uniform" soft thresholding as 

0„(w, s; a) = sign(w n ) max(|u> n | - dy/s^, 0). (23) 

In SUREshrink and PUREshrink [37]-[39] for Gaussian (resp. Poisson) noise reduction, a is 
set to the value that minimizes the corresponding unbiased risk estimate. Similarly, we may select 
a to yield the minimum CURE value according to (19) on the basis of observed data y, resulting 
in a CUREshrink denoising procedure. To comply with the requirements of Theorem 2, we use 
a continuously differentiable approximation to soft thresholding. Figure 4 shows the empirical 
accuracy of CURE as a practical criterion for choosing the best value of a; we have also observed 
a pointwise LET approach, as in (15), to yield comparable denoising results. 

C. Joint inter-Zintra-scale CURE-LET 

To decrease the usual ringing artifacts inherent to orthogonal transform-domain thresholding, 
more sophisticated denoising functions must be considered. In particular, the integration of inter- 
scale dependencies between wavelet coefficients (the so-called "parent-child" relationship) has 
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0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

Value of parameter a Value of parameter a 

Fig. 4. Minimum CURE vs. MMSE threshold selection for the CUREshrink adapted soft thresholding of (23). Left: 
N = 128 x 128 samples, s ~ xJfc), input SNR = 15 dB. Right: N = 256 x 256, s ~ XieW. input SNR = 10 dB. 

already been shown to significantly increase the denoising quality in both AWGN reduction [26] , 
[40]-[42] and Poisson intensity estimation [38]. 

To this end, Fig. 5 shows an example of a group-delay compensated parent p and its child 
w for a particular subband at the first scale of a 2D unnormalized Haar wavelet transform. For 




Fig. 5. A child (left) and its group-delay compensated parent (right) in a particular subband at the first scale of 
a 2D unnormalized Haar wavelet transform. In this 16-color map, the background yellowish value is zero, with the 
most significant coefficients appearing either in blue (negative) or in red (positive). 

the unnormalized Haar wavelet transform, the group-delay compensated parent p = [p n ]i<n<N 
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is simply given by p n = s n+ \ — s ra _i [38]. From Fig. 5, we observe that both the signs and 
the locations of the significant (i.e., large-magnitude) coefficients persist across scale. To take 
advantage of this persistence, we propose the following LET approach, inspired by (15): 

n (w,s;a) = ^flfcinax (l - X k " jw n + ^a fc+2 max (l - X k " ,0)w n 



2 



+ y^ a fe+4max (l - \k ^Z n } S } ,0)pn + VV +6 max (l - A fc 4 ^"/ S \ o)p n . (24) 

^ V 7n(w) / ^ V 7n(P) ' 

Here the function 7 n (u) = \j\f 7 hx Y^ k \u k \e~^ n ~ k ^ 2 / 2 implements a normalized Gaussian smooth- 
ing of the magnitude of its argument; this local filtering accounts for similarities between 
neighboring wavelet, scaling, and parent coefficients. 

The proposed denoising function of (24) thus integrates both the inter- and intra-scale depen- 
dencies that naturally arise in the Haar wavelet transform. It involves a set a of eight parameters 
that can be optimized via least squares, as well as parameters Ai and A2 that can be fixed 
in advance without noticeable loss in denoising quality; we use Ai = 1 and A2 = 9 in all 
experiments below. Considering the processing of all wavelet coefficients in a given subband 
j, (24) reads as 0(w, s; a) = Ylt=i a fc^fc(w, s). The optimal (in the minimum CURE sense) set 
of linear parameters is then the solution of the linear system of equation Ma = c, where 

c = [w T fe (w,s)-4^s-^-l^ di0 fc (w,s) + 8s T d 1 2 2 6> fe (w,s)+ 

4w T {d 2 u e k (w, s) + 3 2 2 2 fc (w, s) - d 2 O k (w, s)) )] i<fc<8 , 
M = [^(w,s) T / (w,s)] 1 < M < 8 . 
Finally, the bias is removed from the lowpass residual subband at scale J as ? J = s J - 2 J K ■ 1. 
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V. Application to Magnitude MR Image Denoising 



In magnitude magnetic resonance imaging, the observed image consists of the magnitudes 
\m n \ of N independent complex measurements m n , where 



U{m n } 



(25) 



Our objective is to estimate the original (unknown) magnitudes |/x n | = y$t{^ n } 2 + 9{A*n} 2 
from their noisy observations \m n \. If we define two iV-dimensional vectors 

x = [|/"n| 2 /0" 2 ]l<n<7V e 



y = [|m„| 2 A^]i<n<iV G R£, 



(26) 



then the data likelihood for y is the product of iV independent noncentral % 2 distributions with 
K = 2 degrees of freedom and noncentrality parameter x n ; i.e., (1) with K = 2. 
We then denoise the magnitude MR image m according to the following steps: 

1) If necessary, estimate the noise variance a 2 using known techniques; 

2) Rescale the squared-magnitude MR image y according to (26); 

3) Apply a CURE-optimized denoising algorithm to obtain an estimate x = f(y) of x; 

4) Fix some A G [0, 1] and produce the final estimate fi of the unknown magnitude MR 
image /x, by applying the following nonlinear rescaling function: 



Vl/„(y)| + (l - A)vW(/ n (y),0) 



(27) 



l<n<iV 

We evaluate denoising performance objectively using three full-reference image quality metrics: 

/VII ii II 2 

1) The standard peak signal-to-noise ratio (PSNR), defined as PSNR = 101og 10 y J^p ; 

2) A contrast-invariant PSNR, defined as CIPSNR = 10 log 10 [r^f^l^jp , where the affine 
parameters a* , b* are given by 



(a* , b*) = argmin \\(afi + b) — /i| 



iVAV - (i T A) 2 ' 

l T /i /i T /i — /i T /i l T /i 

iV/i T /i - (i T A) 2 
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3) The mean of the structural similarity index map (SSIM), a popular visual quality metric 
introduced in [43] (see https://ece.uwaterloo.ca/~z70wang/research/ssim/). 

In all simulated experiments, we have assumed that the variance a 2 of the complex Gaussian 
noise is known. In practice, a reliable estimate can be obtained in signal-free or constant regions 
of the image by moment matching [12], [21] or maximum likelihood techniques [6], [8], [44]. 
When no background is available, more sophisticated approaches can be considered [45]. 

To simulate various input noise levels, several values for a have been selected in the range 
a G [5, 100]. The set of high-quality magnitude MR test images used is shown in Fig. 6, and 
may be obtained from http://bigwww.epfl.ch/luisier/MRIdenoising/TestImages.zip. 



Image 1 (256 X 256) Image 2 (512 X 512) Image 3 (256 X 256) Image 4 (256 X 256) 




Fig. 6. Test set of high-quality magnitude MR images used in the experiments of Section V. 



Image 1 Image 2 




Input PSNR [dB] Input PSNR [dB] 

Fig. 7. PSNR comparisons among pointwise undecimated Haar CURE-LET (Section III-B) and joint intra-/inter-scale 
unnormalized Haar CURE-LET (Section IV-C), shown relative to unnormalized Haar CUREshrink (Section IV-B). 
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Image 1 Image 3 




10 15 20 25 30 10 15 20 25 30 



Input PSNR [dB] Input PSNR [dB] 

Fig. 8. PSNR improvements brought by cycle-spinning the unnormalized Haar wavelet transform. 



Image 1 Image 2 




Fig. 9. Sensitivity to the choice of A in (27); ordinate shows PSNR variations relative to the reference case A = 0. 

A. Comparisons between variants of the proposed approach 

Before comparing our approach with four state-of-the-art MR image denoising methods, we 
first evaluate the performance of the various variants of our approach. In Fig. 7, we compare the 
results obtained using the pointwise CURE-LET thresholding of (15), applied in the undecimated 
Haar wavelet transform domain, and the joint intra-/inter- scale LET denoising function of (24), 
applied to the unnormalized Haar wavelet transform coefficients. These are shown relative to a 
baseline provided by the unnormalized Haar CUREshrink approach of (23). As expected, the 
joint intra-/inter-scale denoising function of (24) outperforms the simple soft-thresholding of (23) 
by 1 — 3dB. Moreover, pointwise thresholding applied in a shift-invariant setting outperforms 
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(by 0.5 — 1.5 dB) a more sophisticated thresholding in a shift- variant one. 

Note that the shift-invariance of the unnormalized Haar wavelet transform can be increased by 
applying the so-called "cycle-spinning" technique [30]. In Fig. 8, we show the PSNR improve- 
ments brought by averaging the results of several cycle-spins (CS). As observed, 16 cycle-spins 
allow a near match in performance to the shift-invariant transform. The cycle-spinning technique 
has the further advantage of being easily implementable in parallel. 

In Fig. 9, we evaluate the sensitivity of the CURE-optimized algorithms with respect to the 
value of parameter A appearing in the final nonlinear reconstruction function of (27). Note that 
this parameter selection is usually not addressed in the literature; the most common choices are 
either A = as in [16] or A = 1 as in [22]. Yet, Fig. 9 indicates that the MMSE choice usually 
lies in between these two extremal values. In all our experiments, we have thus used A = 0.5. 

B. Comparisons to state-of-the-art MR image denoising methods 

As benchmarks for evaluating our CURE-LET approach, we have retained four state-of-the- 
art MR image denoising techniques: two wavelet-based algorithms [21], [22] (code at http: 
//telin.ugent.be/~sanja/), a spatially adaptive linear MMSE filter [12], and an unbiased nonlocal 
means filter specifically designed for MR data [16] (code at http://personales.upv.es/jmanjon/ 
denoising/nlm2d.htm). For each of these methods, we have used the tuning parameters suggested 
in their respective publications and software, except for the linear MMSE filter, where we have 
hand-optimized (in the MMSE sense) the size of the filter support. 

We have considered three CURE-LET variants: 16 cycle-spins of the joint intra-/inter-scale 
thresholding of (24) applied in the unnormalized Haar wavelet transform; the pointwise thresh- 
olding of (15) applied in the undecimated Haar wavelet transform; and the same pointwise 
thresholding applied in a mixed-basis overcomplete transform (an undecimated Haar wavelet 
transform plus an 8 x 8 overlapping block discrete cosine transform, BDCT). Note that a LET 
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spanning several overcomplete bases was previously considered in [34]. 

In Table I, the PSNR results of the various methods are displayed, while the contrast-invariant 
CIPSNR results are reported in Table II. As observed, the proposed CURE-optimized pointwise 
LET applied in a mixed-basis overcomplete transform consistently achieves the best results. The 
average denoising gains are about +2 dB relative to the method of [21], +3 dB compared to [22], 
+3 dB compared to [12], and +0.4 dB compared to [16]. Overall, these CURE-LET methods 
compare very favorably to state-of-the-art approaches, especially under very noisy conditions, 
in which case the signal-dependent nature of the noise is more pronounced. Note that further 
denoising gains are likely to be obtained by considering more sophisticated thresholding rules 
and image-adaptive overcomplete dictionaries (e.g., [46]). 



TABLE I 
PSNR Comparisons 



C 


5 


10 


20 


30 


50 


100 


5 


10 


20 


30 


50 


100 


Image 


Image 1 256 X 256 


Image 2 512 X 512 


Input PSNR 


34.35 


28.09 


21.69 


17.97 


13.28 


6.74 


32.62 


26.48 


20.36 


16.80 


12.32 


6.14 


r pi] 


37.95 


33.46 


29.09 


26.51 


23.20 


18.35 


38.72 


34.13 


29.41 


26.55 


22.85 


17.83 


[22] 


33.13 


30.47 


28.56 


25.13 


18.56 


10.48 


33.79 


32.17 


26.64 


22.00 


16.55 


9.47 


[12] 


37.29 


32.86 


28.10 


25.22 


21.69 


17.33 


37.28 


32.19 


27.31 


24.58 


21.37 


17.64 


[16] 


39.15 


34.97 


30.65 


28.15 


24.79 


19.04 


39.85 


35.98 


31.72 


28.89 


25.08 


19.61 


CURE-LET 
(Haar, CS=16) 


39.00 


34.84 


30.71 


28.15 


24.87 


20.32 


40.02 


35.88 


31.60 


28.95 


25.50 


20.86 


UWT Haar 
CURE-LET 


38.32 


34.48 


30.62 


28.29 


25.23 


20.73 


39.42 


35.95 


32.07 


29.53 


26.27 


21.56 


UWT/BDCT 
CURE-LET 


39.19 


35.00 


30.89 


28.50 


25.38 


21.02 


40.06 


36.14 


32.24 


29.75 


26.37 


21.60 


Image 


Image 3 256 X 256 


Image 4 256 X 256 


Input PSNR 


34.00 


27.77 


21.48 


17.84 


13.32 


7.07 


33.84 


27.73 


21.64 


18.13 


13.71 


7.36 


[21] 


35.88 


31.72 


27.64 


25.21 


22.10 


17.84 


35.75 


31.50 


27.36 


24.99 


22.07 


18.08 


[22] 


29.35 


27.62 


25.89 


24.73 


19.45 


11.65 


28.20 


26.32 


24.63 


24.15 


20.98 


12.31 


[12] 


35.84 


31.24 


26.58 


24.07 


21.00 


17.47 


35.86 


31.19 


26.70 


24.14 


21.06 


17.40 


[16] 


36.25 


32.51 


28.92 


26.67 


23.53 


19.17 


36.16 


32.24 


28.55 


26.43 


23.43 


18.68 


CURE-LET 
(Haar, CS=16) 


36.55 


32.54 


28.72 


26.47 


23.57 


19.67 


36.36 


32.17 


28.23 


25.98 


23.19 


19.41 


UWT Haar 
CURE-LET 


36.43 


32.71 


29.07 


26.88 


23.99 


20.03 


36.20 


32.24 


28.58 


26.49 


23.77 


20.00 


UWT/BDCT 
CURE-LET 


36.88 


33.02 


29.24 


27.01 


24.10 


20.09 


36.56 


32.51 


28.72 


26.59 


23.89 


20.12 



Output PSNRs have been averaged over 10 noise realizations. 
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TABLE II 
CIPSNR Comparisons 



<7 


5 


10 


20 


30 


50 


100 


5 


10 


20 


30 


50 


100 


Image 


Image 1 256 X 256 


Image 2 512 X 512 


input CLraiNK 


34.52 


28.58 


22.72 


19.69 


16.85 


15.16 


34.71 


28.66 


22.59 


19.29 


16.08 


14.13 


[21] 


38.00 


33.57 


29.41 


27.05 


24.02 


19.46 


39.35 


35.15 


31.02 


28.50 


25.13 


20.27 


[22J 


33.55 


30.89 


28.64 


26.29 


22.89 


18.72 


36.41 


33.36 


28.63 


25.88 


22.81 


18.94 


rm 
L12J 


37.30 


32.91 


28.33 


25.57 


22.17 


17.91 


37.81 


33.18 


28.57 


26.06 


23.01 


19.29 


LloJ 


39.15 


34.99 


30.66 


28.16 


24.80 


19.21 


39.93 


36.11 


31.94 


29.17 


25.42 


19.92 


CURE-LET 
(Haar, CS=16) 


39.00 


34.85 


30.82 


28.40 


25.30 


20.75 


40.29 


36.37 


32.50 


30.13 


26.98 


22.46 


UWT Haar 
CURE-LET 


38.32 


34.48 


30.64 


28.36 


25.41 


21.03 


39.56 


36.17 


32.45 


30.07 


27.08 


22.71 


UWT/BDCT 
CURE-LET 


39.20 


35.01 


30.91 


28.58 


25.60 


21.36 


40.30 


36.45 


32.68 


30.34 


27.22 


22.79 


Image 


Image 3 256 X 256 


Image 4 256 X 256 


Input CIPSNR 


34.33 


28.31 


22.24 


18.88 


15.44 


13.06 


34.07 


28.01 


22.13 


19.06 


16.07 


14.07 


[21] 


36.01 


31.88 


28.00 


25.76 


22.90 


18.83 


35.94 


31.77 


27.79 


25.51 


22.64 


18.53 


[22] 


30.25 


28.49 


26.35 


24.98 


22.02 


18.13 


28.96 


26.93 


25.07 


24.46 


22.34 


18.25 


[12] 


35.86 


31.38 


26.88 


24.58 


21.65 


18.17 


35.93 


31.35 


26.93 


24.41 


21.34 


17.68 


[16] 


36.27 


32.51 


28.93 


26.67 


23.56 


19.21 


36.17 


32.26 


28.57 


26.46 


23.48 


18.72 


CURE-LET 
(Haar, CS=16) 


36.58 


32.59 


28.86 


26.73 


24.03 


20.29 


36.42 


32.27 


28.44 


26.27 


23.55 


19.67 


UWT Haar 
CURE-LET 


36.48 


32.74 


29.12 


26.97 


24.19 


20.44 


36.25 


32.28 


28.66 


26.61 


23.96 


20.24 


UWT/BDCT 
CURE-LET 


36.93 


33.07 


29.29 


27.11 


24.31 


20.55 


36.61 


32.56 


28.82 


26.73 


24.09 


20.42 



Output CIPSNRs have been averaged over 10 noise realizations. 



Fig. 10 presents a visual comparison of the various algorithms considered. As observed, the 
two CURE-LET denoising results offer a good balance between noise suppression, reasonable 
denoising artifacts, and fine structure preservation. This subjective observation is confirmed by 
the higher SSIM scores obtained by the proposed denoising approach. 

In Table III, we report the computation time of the various algorithms considered. All have 
been executed on Matlab R2010a running under Mac OS X equipped with a 2.66GHz Intel 
Core 2 Duo processor. As observed, the proposed CURE-LET algorithms are quite fast, taking 
1-10 s to denoise a 256 x 256 image. When applied within an undecimated filterbank transform, 
most of the computational load is dedicated to the independent reconstructions of the processed 
subbands and their corresponding first and second order derivatives (see (16)). 
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Fig. 10. Visual Comparisons. (A) Zoom in image 1. (B) Noisy version: SSIM = 0.319. (C) Denoised by [21]: 
SSIM = 0.734. (D) Denoised by [22]: SSIM = 0.726. (E) Denoised by [12]: SSIM = 0.653. (F) Denoised 
by [16]: SSIM = 0.719. (G) Denoised by Haar CURE-LET (CS=16): SSIM = 0.800. (H) Denoised by UWT/BDCT 
CURE-LET: SSIM = 0.796. 

C. Denoising of a magnitude MR knee image 

We have also applied our CURE-LET denoising algorithms to an actual magnitude MR image 
of the knee. This 512x512 16-bit raw image has been acquired on a Siemens 1.5 Tesla Magnetom 
Sonata MR system, following a sagittal T2-weighted protocol. The standard deviation of the 
complex Gaussian noise has been estimated from a signal-free region S of the squared data, as 
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TABLE III 

Computation times for MR image denoising techniques. 



Method 


Image size 


Method 


Image size 


256 x 256 


512 X 512 


256 x 256 


512 x 512 


[21] 1 


0.2 


0.9 


Haar CURE-LET (CS=1) 


0.3 


0.8 


[22] 


0.6 


2.5 


Haar CURE-LET (CS=16) 


3.7 


12.1 


[12] 


0.1 


0.5 


UWT CURE-LET 


1.6 


9.6 


[16] 


55.8 


248.4 


UWT/BDCT CURE-LET 


9.7 


43.9 



Computation times have been averaged over 10 runs; [12], [16] do not use pre-compiled MEX files. 



a = y ^ J2neS l m "l 2 > an ^ subsequently treated as known. 

Fig. 1 1 shows the denoising results of the various CURE-LET algorithms. As observed, the 
noise is efficiently attenuated and the contrast is significantly improved, owing to a proper 
reduction of the signal-dependent bias introduced by the noise. 

VI. CONCLUSION 

In this article we have derived a noncentral chi-square unbiased risk estimate (CURE), and 
applied it to the problem of magnitude MR image denoising, where the squared value of each 
pixel comprises an independent noncentral chi-square variate on two degrees of freedom. Our 
approach can be used to optimize the parameters of essentially any continuously differentiable 
estimator for this class of problems, and here we have focused our attention on transform-domain 
algorithms in particular. 

In this vein, we first developed a pointwise linear expansion of thresholds (LET) estimator 
applied to the coefficients of an arbitrary undecimated filterbank transform. We then considered 
the specific case of the unnormalized Haar wavelet transform, a multiscale orthogonal transform 
allowing for the derivation of subband-dependent CURE denoising strategies. We also introduced 
a subband-adaptive joint inter-/intra-scale LET that outperforms a simpler estimator similar to 
soft thresholding. 
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(A) (B) 




Fig. 11. Denoising of a magnitude MR image of the knee. (A) Raw 16-bit data: a = 67.61. (B) Denoised via Haar 
CURE-LET (16 cycle spins). (C) Via UWT Haar CURE-LET. (D) Via UWT/BDCT CURE-LET. 

We then applied our proposed CURE-optimized algorithms to test images artificially degraded 
by noise, and observed them to compare favorably with state-of-the-art techniques, both quan- 
titatively and qualitatively. Finally, we showed an example of denoising results obtained on an 
actual magnitude MR image, in order to show the practical efficacy of our approach to MR 
image denoising via chi-square unbiased risk estimation. 
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